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ABSTRACT 


INTERACTIVE REAL TIME FLOW 
SIMULATIONS 

Ideen Sadrehaghighi 
Surendra N. Tiwari 

An interactive real time flow simulation technique is developed for an un- 
steady channel flow. A finite-volume algorithm in conjunction with a Runge-Kutta 
time stepping scheme has been developed for two-dimensional Euler equations. A 
global time step has been used to accelerate convergence of steady-state calculations. 
A raster image generation routine has been developed for high speed image transmis- 
sion which allows user to have direct interaction with the solution development. In 
addition to theory and results, the hardware and software requirements are discussed. 
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= vector of conservative variables 
= flux vectors for coordinate directions 
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= angle of attack 
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Chapter 1 

INTRODUCTION 


Since an interactive design process is one of the ultimate goals of Computa- 
tional Fluid Dynamics (CFD), real time flow simulation with direct user interaction 
is an ideal approach in achieving this goal [1]*. This process requires a supercomputer 
with vast amount of memory, an extremely high bandwidth communication network 
and highly capable graphic workstations. Even with today’s rapid advances in the 
supercomputer developments, some of the components are not fast or large enough 
for a realistic three-dimensional problem. However, real time simulation of medium 
size two-dimensional flow problems (Euler equations) are possible today. 

Accurate simulations are critical to the development and understanding 
of highly unsteady flow. For the reasons outlined above, a simple unsteady two- 
dimensional channel type flow has been studied here. The relative motion of the 
shocks and other strong gradients have been examined as the solution being, com- 
puted. Once the solution has been examined, the user will inform the flow solver 
to take different action or to continue on the existing course. This requires large 
amount of data to be examined by the user, and small amount of information in 
form of instruction, to be send back to the mainframe computer. Consequently, this 
procedure requires a fast mainframe computer, an extremely fast network for data 
communication and a relatively fast workstation. 

The computational process is initiated on the mainframe computer under 

* The numbers in brackets indicate references. 
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interactive control by the user at a workstation. The equations of motion are in- 
tegrated step by step on the mainframe computer. The choice of a variable to be 
visualized is instructed from the workstation and a separate rasterization program 
computes a raster image of the variable. The image is transmitted over the communi- 
cation link to a raster display device to be viewed by the user. Images are continuously 
created, transmitted, and if desired, stored on a recording device. Having viewed the 
images on the raster display device, the user responds. If mainframe’s computational 
and communication rates are sufficient, a real time interaction can be achieved. This, 
of course , depends on the size of the problem under investigation (i.e. number of 
grid points, equations of motion, and solution technique) . 


Chapter 2 

GRID GENERATION 

A suitable transformation for a channel type flow can be obtained using an 
analytical function. A simple case would be to use uniform spacing in x direction. 
For y direction, the spacing is obtained by [2], 

(g + 2y)( (/? + 1 )/«? - i - 0 + 2, 

y ~ (2>r + l)(i + [(/J + 1)/(/9 - )) 

where h is the height of channel and ij is the parameter which controls the grid 
clustering in y direction. Variable y corresponds to normal coordinate in the com- 
putational domain. The stretching parameter, /?, is related approximately to the 
non-dimensional boundary layer thickness (S/h) by 

)J = (l-£)* , 0 < I < 1 (2.2) 

The amount of stretching for various values of 6/h is illustrated in Fig. 2.1. 
For this transformation, if rj = 0 the mesh will be refined near y = h , whereas, if 
T) = 0.5 the mesh \yill be refined equally near y = 0 and y = h . For this study, values 
of rj = 0.5 and /?= 0.30 been chosen and the resulting grid is shown in Fig. 2.2. 
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Fig. 2.2 Computational grid for a channel flow 







Chapter 3 


EQUATIONS OF MOTION AND 
METHOD OF SOLUTION 

3.1 Governing Equations 

The governing equations are time-dependent Euler equations [3], The un- 
steady two-dimensional equations for a compressible perfect gas can be written in an 
integral form for a region fl with boundary as 

~ J J n Q dxdy + JjFdy - Gdx) = 0 (3.1) 

where vectors Q , F , and G, are 

■ p 1 pu ] r pv 

Q _ P U F= P^+P Q_ pUV 

pv puv pv 2 + p 

. E . . {E -f p)u . . (E + p)v 

The total energy, E, is defined using the ideal gas relation as 

E = ( 7 = 7 ) + W + < 3 - 3 ) 

where 7 is the ratio of specific heats. 
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3.2 Finite- Volume Scheme 

The governing equations in the integral form are applied to an arbitrary 
quadrilateral and the line integrals are approximated with the midpoint rule. The 
discretized result is 


J t (Si, Qi.>) + = 0 (3.4) 

where L is defined as the spatial discretization operator and S,j is area of the cell. 
The components of Q , j represent the cell-averaged quantities and obtained by the 
mean values of the fluxes crossing the cell. 

In order to overcome the instabilities associated with central differencing, a 
fourth order dissipative model is used. In this study, the dissipative model developed 
by Jameson, Schmidt, and Turkel [4] has been used. The finite-volume formulation, 
Eq. (3.4) is now expressed as 

Qij) + £Qi,j - DQij = 0 (3.5) 

where D represents the artificial dissipation operator. This operator can be written 
in the following way 

DQtj = D x Qij + D y Q tJ (3.6) 

where D x Qij and D y Qij are the contributions from each of the coordinate directions. 
In conservation form 

D x Qi,j = d, +1 / 2 j — d«-i/2,j (37) 

A/Qui = ^iJ+i/2 ~ d»,i- 1/2 (3.8) 

where the terms on the right hand side have the form 
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S' 

d;+i/2,j = — '[ £ !+ ) |/2,j(Qi+ij - Q.’.j) — e !+\/2j(Q'+2j — 3Qi + i,j + 3Q tJ - Q;_ ltJ )] 

(3.9) 

The coefficients and are determined from the pressures gradients as 


f . | Pj+ij - 2Pjj ± Pj-ij | 

P i+lJ + 2 Pij + P t _ u 

e !+ 1/2,2 = <* { 2 ) max(v i+u v iti ) 

e !+ i/2,i = max (°> (« (4) ~ e\+x/2j)) 

where typical values of the constants a< 2) and a (4) are ± and respectively. 


(3.10) 

(3.11) 

(3.12) 


3.3 Time- Stepping Scheme 

A modified four stage Runge-Kutta time stepping scheme was used to ad- 
vance the solution in time. Due to its explicit character, this scheme is very simple 
and flexible . Hence, it is ideal for interactive computation purposes. At time level n, 


Q (0) = Q (n) (3.13) 

Q (1) = Q (0) -a,A<rtQ<°> (3.14) 

Q (2) = Q(°) - Q2 AtRQM (3.15) 

Q (3) = Q (0) - a 3 AtRQ^ (3.16) 

QD) _ Q(°) _ AtRQO) (3.17) 

Q("+l) = q«> (3. 18 ) 

where on the (q+l)st stage 


RQM = ^(IQW - DQ(°)) 


(3.19) 
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and 


111 

«i = ^ » ° 2 = 3 ’ ° 3 = 2 (3.20) 

3.4 Boundary Conditions 

The boundary conditions are implemented by applying a slip velocity con- 
dition at the solid walls [5). Also, the stagnation enthalpy is assumed to be constant 
along the solid surface and equal to that of the free-stream, i.e., 

ht = h + -(u 2 + u 2 ) = h too (3.21) 

For solid wall, (h t )„ in u, the free-stream values(u=l and v= 0 ) arc used, i.e., 

{ht) W aii = J^F.s. + 2 (3.22) 

Because of the supersonic nature of the flow, the outflow boundary is determined by 
extrapolating from interior points. 



Chapter 4 


SOFTWARE AND HARDWARE 
CONSIDERATIONS 


The implemented software and hardware are discussed in [1], Currently, 
the mainframe computer is a CRAY-2 which generally performs at 200-250 M FLOPS 
(Million Floating Point Operations Per Second) . The mainframe is connected to 
UltraNet Graphics Display Device (UGDD), which is a network based frame buffer 
through a High-Speed Channel (HSC) . The channel can support transfer rates up 
to 100 MBytes/sec. The frame buffer is a high-speed display system supporting high 
resolution monitors. The UGDD contains two memory arrays. While one array is 
displayed, the other is updated over the network. The user determines which array 
to display and controls the situation over the network. The UGDD supports a single 
user at a time and may be used to display color images at animation rates. The 
frame buffer can display RGB (Red, Green, Blue) images up to 1280 pixels horizon- 
tally by 1024 Scan-lines vertically. Figure 4.1 shows a pictorial representation of the 
hardware setup for a single hub (UltraNet 1000). Figure 4.2 shows the configuration 
using a CRAY supercomputer in conjunction with several mini-supercomputers and 
workstations. 

There are several pieces of software required for a real time CFD simulation. 
They are : a grid generator, a flow solver, an image rasterizer, and an image manip- 
ulator. Each software component has been written and applied separately. However, 
for a true interactive process, these softwares must be integrated into one working 


9 


10 


unit or be presented in parallel. 

With the current CFD technology, it is possible to solve the Euler equations 
at 1 67 — — r [61. For a grid of 70 X 30, it will take one second for 62 
iterations. This will result in 62 frames per second. 

A rasterization program called PLOTD is developed to produce a raster 
image from the solution. A complete listing of the source code is provided in Ap- 
pendix A. This software converts a specified variable defined 011 a surface to a raster 
RGB images (1024 X 768). The required CPU time for this step depends on the 
size of the grid and the complexity of the image. This step may take from few mi- 
croseconds to a couple of seconds. Consequently, the image must not be very complex. 













Chapter 5 

RESULTS AND DISCUSSION 


The unsteady two-dimensional Euler equations are solved for a channel flow. 
In order to test the flow solver, the solution has been obtained by keeping all the flow 
parameters constant. The results are satisfactory and the shocks are captured very 
accurately. Figures 5.1 through 5.4 show pressure and Mach number contours at 
different iterations for different angles of attack. A more comprehensive test can be 
achieved by changing one of the flow conditions as the solution evolves. An ideal and 
relatively simple case would be rotating the angle of attack. The angle of attack can 
be expressed as 


a = a min -f (a mor - Q min )stn(u;ir) (5.1) 

where a mm and ct max are intial and final angles of attack, and u> is the specified 
frequency. Initially, the solution starts with a uniform flow everywhere (impulse 
start). At this initial stage of the computation, none of the flow features lias been 
altered. This, in fact, is the continuation of the case described earlier for testing the 
flow solver. After some iterations, once the flow features are established, the angle 
of attack is allowed to rotate using Eq. (5.1). It takes between 3-5 cycles for flow 
properties to reach a cyclic behavior. 

Each frame takes at least 0.5 second of CPU time, which is relatively long 
time for a multiuser machine. One way to alleviate this is to take advantage of the 
parallel capability of the CRAY 2. This can be accomplished by allowing the flow 


12 



13 


solution, rasterization and image transmission to be performed on different processors. 





























after 95 iterations 



after 45 iterations 


Fie- 5.4 Mach number contours for different iterations (o = 15*0 , = 2.0) 


REFERENCES 


1. Abolhassani, J. S. and Everton, E. L., “ Interactive Grid Adaption , ” The Third 
International Congress of Fluid Mechanics , Vol. I, January 1990, pp. 309-317. 

2. Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., Computational Fluid 
Mechanics and Heat Transfer . Hemisphere Publishing Corporation, New York, 
1984. 

3. Lavantc, E. V., Hruns, R. L., Sanetrik, M. D., and Lam, T., “ Numerical Analysis 
of Flow About a Total Temperature Sensor ,” AIAA Journal, March 1989, pp 692 
- 714. 

4. Jamson, A., Schmidt, W. and Turkel, E., “ Numerical Solutions of the Eu- 
ler Equations By Finite- Volume Using Runge-Kutta Time-Stepping Scham.es ,” 
AIAA Paper 81-1259, June 1981. 

5. IIofFmann, K. A., Computational Fluid Dynamics For Engineers , Engineering 
Education System, Austin, TX ,1989. 

6. Rumsey, C. L. and Anderson, W. K., u Some Numerical and Physical Aspects of 
Unsteady Navier-Stokes Computations Over Airfoils Using Dynamic Meshes," 
AIAA Paper 88-0329, January 1988. 


18 


APPENDIX A 

RASTER IMAGE GENERATION ROUTINE : PLOTD 


Following subroutines are developed for a raster image generation. Main 
arguments are defined as: 


*,y 

d 

n 

m 

is,ie 

js,je 

nnc 


= grid coordinates 
= flow variable to be contoured 
= number of points in x direction 
= number of points in y direction 

= desired start and end points for contouring in I-direction 
= desired start and end points for contouring in J-direction 
= number of contour levels 
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subroutine conpl t ( x,y,d,m,n,is,ie, js, je,a,b, nnc ) 
c 

PARAMETER (NPMAX= 10000 , NPSMAX=100 , NSMAX=100 ) 
c 

dimension x(m, n) # y(m, n) , d(m, n) 
c 

common /extrem/ xmin, xmax, ymin, ymax, dmin, 

real xmin, xmax, ymin, ymax, dmin, 

c 

common /maping/ scalex, ofsetx, scaley, ofsety, scaled, 

real scalex, ofsetx, scaley, ofsety, scaled, 

c 

dimension ixl ( npsmax , 3 ) , iyl ( npsmax , 3 ) 
dimension ix2 ( npsmax , 3 ) , iy2 ( npsmax , 3 ) 
c 

logical ind( npsmax, 2 ) , flag (npsmax) 
c 

c contour plot 

c written by: Eric L. Everton & Jamshid S. Abolhassani 

c 

nc=abs ( nnc ) 
c 

c check for zero range 

c 

if ( dmax . ne . dmin ) then 
c 

c compute delta 

c 

delf=(dmax-dmin)/f loat (nc ) 
c 

else 

c 

c zero range, issue error message and stop 

c 

write(*,*) 'maximum and minimum are equal' 
stop 
c 

endi f 
c 

c set starting contour level 

c 

cl=dmin 

c 

c set error flag to false 

c 

do 100 i=is,iel 
c 

f lag( i ) = . false . 
c 

100 continue 
c 

c plot contours 

c 

150 continue 
c 

c checks for contour level out of range 


dmax 

dmax 

of setd 
of setd 
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c 

if ( cl . It . dmin ) go to 150 
c 

if ( cl . gt . dmax ) return 
c 

c convert contour level to look up table (lut) index 

c 

icolor=cl*scaled+of setd 
c 

c checks for lut index out of range 

c 

if(icolor.lt.int(a)) icolor=int(a) 
c 

if (icolor . gt . int(b) ) icolor=int(b) 
c 

c set current lut index (color) 

c 

call ufbcol ( icolor ) 
c 

c set do loop end controls 

c 

iel=ie- 1 
jel=je-l 
c 

c upper triangle do loop 

c 

do 200 j=js,jel 
c 

c reset line plot flags 

c 

do 250 k=l , 2 
do 250 i=is , iel 
c 

ind(i,k)=. false, 
c 

250 continue 
c 

c vectorized do loop 

c 

do 300 i=is, iel 
c 

c initialize temporary variables 

c 

xl=x( i , j ) 

yi=y(i ,j ) 

f l=d( i ,j ) 
x2=x(i+l,j ) 
y2=y ( i+1 , j ) 
f2=d( i+1 , j ) 
x3=x( i+1 , j+1 ) 
y3=y( i+1 , j+1 ) 
f3=d(i+l, j+1) 
c 

c check to see if contour line intersects this triangle 

c 

if ( (cl . le .max( f 1 , f2 , f3 ) ) . and . (cl . ge . min( f 1 , f2 , f 3 ) ) ) then 
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c 

c check to see if contour line intersects point 1 

c 

if (cl.eq.fi) then 
c 

c check to see if contour line intersects point 2 

c 

if (cl.eq.f2) then 
c 

c check to see if contour line intersects point 3 

c 

if (cl.eq.f3) then 
c 

c contour line intersects all three verticies 

c 

ixl ( i , 1 )=of setx+scalex*xl 
iyl ( i , l)=ofsety+scaley*yl 
ix2 ( i , 1 )=of setx+scalex*x2 
iy2 ( i , 1 ) =of sety+scaley*y2 
ind( i , 1 )= . true . 
ixl(i,2)=ofsetx+scalex*x2 
iyl ( i , 2 )=of sety+scaley*y2 
ix2 ( i , 2 )=of setx+scalex*x3 
iy2 ( i , 2 )=of sety+scaley*y3 
ind( i , 2 )= . true . 
ixl(i,3)=ofsetx+scalex*x3 
iyl ( i , 3 )=of sety+scaley*y3 
ix2(i,3)=ofsetx+scalex*xl 
iy2(i,3)=ofsety+scaley*yl 
c 

else 

c 

c contour line intersects point 1 and point 2 

c 

ixl ( i , 1 )=of setx+scalex*xl 
iyl(i, l)=ofsety+scaley*yl 
ix2(i,l)=ofsetx+scalex*x2 
iy2 ( i , 1 )=of sety+scaley*y2 
ind( i , 1 )= . true . 
c 

endi f 
c 

c contour line did not intersect point 2 

c check to see if contour line Intersects point 3 

c 

else if (cl.eq.f3) then 
c 

c contour line intersects point 1 and point 3 

c 

ixl ( i , 1 )=of setx+scalex*xl 
iyl(i, l)=ofsety+scaley*yl 
ix2(i,l)=ofsetx+scalex*x3 
iy2 ( i , 1 )=of sety+scaley*y3 
ind( i , 1 )=. true . 


c 


else 



check to see if contour line intersects edge 
between point 2 and point 3 

if ( (cl-f2)*(cl-f3) . le.O. ) then 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


contour line intersects point 1 and 
edge between point 2 and point 3 


ixl ( i , l)=ofsetx+scalex*xl 
iyl ( i , l)=ofsety+scaley*yl 

ix2(i, l)=ofsetx+scalex*(x2+(x3-x2)*(cl-f2)/(f3-f2) ) 
iy2(i,l)=ofsety+scaley*(y2+(y3-y2)*(cl-f2)/(f3-f2) ) 
ind( i , 1 ) = . true . 


endif 


endif 


contour line did not intersect point 1 

check to see if contour line intersects point 2 

else if (cl.eq.f2) then 

contour line intersects point 2 

ixl(i,l)=ofsetx+scalex*x2 
iyl( i , 1 )=of sety+scaley*y2 

check to see if contour line intersects point 3 
if (cl.eq.f3) then 

contour line intersects point 2 and point 3 

ix2 ( i , 1 )=of setx+scalex*x3 
iy2 ( i , 1 )=of sety+scaley*y3 
ind ( i , 1 )= . true . 

else 

check to see if contour line intersects edge 
between point 1 and point 3 

if ( (cl-fl)*(cl-f3) . le.O. ) then 

contour line intersects point 2 and 
edge between point 1 and point 3 

ix2 ( i , 1 )=of setx+scalex* (xl+ (x3-xl )*(cl-fl)/(f3-fl) ) 
iy2 ( i , 1 )=of sety+scaley* (yl+(y3-yl)*(cl-fl)/(f3-fl) ) 
ind( i , 1 )= . true . 

endif 


c 


endif 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


check to see if contour line intersects point 3 

else if (cl.eq.f3) then 

check to see if contour line intersects edge 
between point 1 and point 2 

if ( (cl-f 1 )*(cl-f2 ) . le.O. ) then 

contour line intersects point 3 and 
edge between point 1 and point 2 

ixl(i,l)=ofsetx +scalex*x3 
iyl ( i , 1 )=of sety+scaley*y3 

ix2(i, 1 )=of setx+ scalex*(xl+(x2-xl)*(cl-fl)/(f2-fl) ) 
iy2 ( i , 1 )=of sety+scaley* (yl+(y2-yl)*(cl-fl)/(f2-fl)) 
ind( i f 1 )= . true . 

endif 

else 

contour line does not intersect any of the verticies 
check to see if contour line intersects edge 
between point 2 and point 3 

if ( (cl-f2)*(cl-f3) . le.O. ) then 

contour line intersects edge 
between point 2 and point 3 

ixl ( i f l)=ofsetx+scalex*(x2+(x3-x2)*(cl-f2)/(f3-f2) ) 
iyl ( i , 1 )=of sety+scaley* ( y2+ ( y3-y2 ) * ( cl-f 2 )/( f 3- f 2 ) ) 

check to see if contour line intersects edge 
between point 1 and point 3 

if ( (cl-fl)*(cl-f3) .le.O. ) then 

contour line intersects edges 
between point 2 and point 3 and 
between point 1 and point 3 

ix2(i, l)=ofsetx+scalex+(xl*(x3-xl)*(cl-fl)/(f3-fl) ) 
iy2 ( i , 1 )=of sety+scaley* (ylMy3-yl)*(cl-fl)/(f3-fl)) 
ind( i , 1 )= . true . 

contour line does not intersect edge 
between point 1 and point 3 

check to see if contour line intersects edge 
between point 1 and point 2 


c 

c 


else if ((cl-fl)*(cl-f2) . le.O. ) then 
contour line intersects edges 



25 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


between point 2 and point 3 and 
between point 1 and point 2 

ix2 ( i , l)=ofsetx+scalex*(xl+(x2-xl)*(cl-fl)/(f2-fl) ) 
iy2 ( i , 1 )=of sety+scaley* (yl+(y2-yl)*(cl-fl)/(f2-fl) ) 
ind( i , 1 ) = « true . 

must be a problem 

else 

set error flag 
f lag( i )= . true . 
endif 

contour line does not intersect edge 
between point 2 and point 3 

check to see if contour line intersects edge 
between point 1 and point 3 

else if ((cl-fl)*(cl-f3) • le.O. ) then 

contour line intersects edge 
between point 1 and point 3 

ixl(i l l)=ofsetx+scalex*(xl+(x3-xl)*(cl-fl)/(f3-fl)) 

iyl(i, l)=ofsety+scaley*(yl+(y3-yl)*(cl-fl)/(f3-fl) ) 

check to see if contour line intersects edge 
between point 1 and point 2 

if ( (cl-fl)*(cl-f2) . le.O. ) then 

contour line intersects edges 
between point 1 and point 3 and 
between point 1 and point 2 

ix2(i,l)=ofsetx+scalex*(xl+(x2-xl)*(cl-fl)/(f2-fl)) 
iy2(i,l)=ofsety+scaley*(yl+(y2-yl)*(cl-fl)/(f2-fl)) 
ind( i , 1 )= . true . 

must be a problem 

else 

set error flag 
f lag( i )= . true . 
endif 


contour line did not intersect edges 
between point 2 and point 3 and 
between point 1 and point 3 
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c check to see if contour line intersects edge 

c between point 1 and point 2 

C else if ( (cl-fl)*(cl-f2) . le.O. ) then 

c # 

c contour line only intersects edge 

c between point 1 and point 2 

c must be a problem 

c set error flag 

c 

f lag( i )=. true . 
c 

endif 


c 

endif 

c 

endif 

c 

300 continue 
c 

do 350 i=is,iel 
c 

c check to see if an error occured 

c 

if ( f lag( i ) ) then 
c 

c issue error message and stop 

c 

write(*,*) ' There is something fishy about conplt!' 
stop 
c 

endi f 
c 

350 continue 
c 

do 360 i=is,iel 
c 

c check to see if line is to be plotted 

c 

if (ind(i,l)) then 
c 

c plot a single line 

c 

call ufblin(ixl(i,l) , iyl(i, 1) . ix2 ( 1 . 1 ) . iy2 ( i , 1 ) ) 
c 

c check for more lines 

c 

if (ind(i,2)) then 
c 

c plot two more lines 

call ufblin(ixl(i, 2 ),iyl(i, 2 ),ix 2 (i, 2 ),iy 2 (i, 2 )) 
call ufblin(ixl(i,3),iyl(i,3).ix2(i,3),iy2(i,3)) 
c 

endif 

c 
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endif 

c 

360 continue 
c 

200 continue 
c 

c lower triangle do loop 

c 

do 400 j=j s , j el 
c 

c reset line plot flags 

c 

do 450 k=l , 2 
do 450 i=is , iel 
c 

ind( i , k) = . false . 
c 

450 continue 
c 

c vectorized do loop 

c 

do 500 i=is, iel 
c 

c initialize temporary variables 

c 

xl=x( i ,j ) 
yl=y(i ,j ) 
f l=d( i ,j ) 
x3=x( i + 1 , j + 1 ) 
y3=y(i+l, j+1) 
f3=d( i+1 , j+1 ) 
x4=x ( i , j + 1 ) 
y4=y(i ,j+l) 
f4=d( i ,j+l) 
c 

c check to see if contour line intersects this triangle 

c 

if ( ( cl . le . max( f 1 , f 3 , f 4) ) . and . ( cl . ge . min( f 1 , f3 , f4 ) ) ) then 
c 

c check to see if contour line intersects point 1 

c 

if (cl.eq.fi) then 
c 

c check to see if contour line intersects point 3 

c 

if (cl.eq.f3) then 
c 

c check to see if contour line Intersects point 4 

c 

if (cl.eq.f4) then 
c 

c contour line intersects all three verticies 

c 

ixl ( i , 1 )=of setx+scalex*xl 
iyl ( i , 1 )=of sety+scaley*yl 
ix2 ( i , 1 )=of setx+scalex*x3 
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iy2 ( i , 1 ) =of sety+scaley*y3 
ind( i , 1 )= . true . 
ixl(i,2)=ofsetx+scalex*x3 
iyl ( i , 2 )=of sety+scaley*y3 
ix2 ( i , 2 )=of setx+scalex*x4 
iy2 ( i , 2 )=of sety+scaley*y4 
ind( i , 2 )= . true . 
ixl ( i , 3 )=of setx+scalex*x4 
iyl ( i , 3 )=of sety+scaley*y4 
1x2 ( i , 3 )=of setx+scalex*xl 
iy2 ( i , 3 )=of sety+scaley*yl 
c 

else 

c 

c contour line intersects point 1 and point 3 

c 

ixl ( i , 1 )=of setx+scalex*xl 
iyl(i, 1) =ofsety+ scale y*yl 
ix2(i, 1 )=of setx+scalex*x3 
iy2(i, l)=ofsety+scaley*y3 
ind( i , 1 )= . true . 
c 

endif 


c 

c contour line did not intersect point 3 

c check to see if contour line intersects point 4 

c 

else if (cl.eq.f4) then 
c 

c contour line intersects point 1 and point 4 

c 

ixl(i, l)=ofsetx+scalex*xl 
iyl ( i , 1 )=of sety+scaley*yl 
ix2 ( i , 1 )=of setx+scalex*x4 
iy2 ( i , 1 )=of sety+scaley*y4 
ind( i , 1 )= . true . 


c 

c 

c 

c 

c 

c 

c 

c 


c 


check to see if contour line intersects edge 
between point 3 and point 4 

if ( (cl-f3)*(cl-f4) . le.O. ) then 

contour line intersects point 1 and 
edge between point 3 and point 4 

ixl(i, l)=ofsetx+Bcalex*xl 
iyl(i, l)=ofsety+scaley*yl 

ix2(i,l)=ofsetx+scalex*(x3+( x4-x3 )*(cl-f3)/(f4-f3) ) 
iy2 ( i , 1 )=of sety+sca ley* ( y3 + ( y4-y3 ) * ( cl-f 3 )/( f 4-f 3 ) ) 
ind( i , 1 )= . true . 

endif 


c 


endif 



contour line did not intersect point 1 

check to see if contour line intersects point 3 

else if (cl.eq.f3) then 
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c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


contour line intersects point 3 

ixl ( i , 1 )=of setx+scalex*x3 
iyl ( i t 1 )=of sety+scaley*y3 

check to see if contour line intersects point 4 
if (cl.eq.f4) then 

contour line intersects point 3 and point 4 

ix2(i, l)=ofsetx+scalex*x4 
iy2 ( i f 1 )=of sety+scaley*y4 
ind( i, 1)=. true . 

else 

check to see if contour line intersects edge 
between point 1 and point 4 

if ((cl-fl)*(cl-f4) .le.O. ) then 

contour line intersects point 3 and 
edge between point 1 and point 4 

ix2 ( i , 1 )=of setx+scalex* ( xl+ ( x4-xl )*(cl-fl)/(f4-fl) ) 
iy2 ( i , 1 )=of sety+scaley* ( yl + ( y4-yl )*(cl-fl)/(f4-fl)) 
ind( i , 1 )=. true . 

endif 

endif 

check to see if contour line intersects point 4 

else if (cl.eq.f4) then 

check to see if contour line intersects edge 
between point 1 and point 3 

if ( (cl-fl)*(cl-f3) . le.O. ) then 

contour line intersects point 4 and 
edge between point 1 and point 3 

ixl ( i , 1 )=of setx+scalex*x4 
iyl ( i t 1 )=of sety+scaley*y4 

ix2 ( i t 1 )=of setx+scalex* (xl+(x3-xl)*(cl-fl)/(f3-fl)) 
iy2 ( i , 1 )=of sety+scaley* (yl+(y3-yl)*(cl-fl)/(f3-fl)) 
ind( i f 1 )=. true . 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


endif 

else 

contour line does not intersect any of the verticies 
check to see if contour line intersects edge 
between point 3 and point 4 

if ( (cl-f3)*(cl-f4) • le.O. ) then 

contour line intersects edge 
between point 3 and point 4 

ixl(i,l)=ofsetx+scalex*(x3+( x4-x3 )*(cl-f3)/(f4-f3) ) 
iyl ( i , 1 )=of sety+ scale y* ( y3+ ( y4-y3 )*(cl-f3)/(f4-f3) ) 

check to see if contour line intersects edge 
between point 1 and point 4 

if ( (cl-fl)*(cl-f4) .le.O. ) then 

contour line intersects edges 
between point 3 and point 4 and 
between point 1 and point 4 

ix2(i > l)=ofsetx+scalex*(xl+( x4-xl )*(cl-fl)/(f4-fl) ) 
iy2 ( i , 1 )=of sety+scaley* ( yl+ ( y4-yl )* (cl-fl)/(f4-fl)) 
ind( i , 1 )=. true . 

contour line does not intersect edge 
between point 1 and point 4 

check to see if contour line intersects edge 
between point 1 and point 3 

else if ( (cl-f 1 )* (cl-f3 ) . le . 0. ) then 

contour line intersects edges 
between point 3 and point 4 and 
between point 1 and point 3 

ix2(i > l)=ofsetx+scalex*(xl+(x3-xl)*(cl-fl)/(f3-fl)) 
iy2(i,l)=ofsety+scaley*(ylMy3-yl)*(cl-fl)/(f3-fl)) 
ind( i , 1 )= . true . 

must be a problem 

else 

set error flag 
f lag( i )= . true . 
endif 


c 

c 


contour line does not intersect edge 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


between point 3 and point 4 

check to see if contour line intersects edge 
between point 1 and point 4 

else if ( (cl-fl)*(cl-f4) . le.O. ) then 

contour line intersects edge 
between point 1 and point 4 

ixl(i,l)=ofsetx+scalex*(xl+( x4-xl )*(cl-fl)/(f4-fl) ) 
iyl ( i , 1 )=of sety+scaley* ( yl + ( y4-yl ) * ( cl-f 1 )/( f 4-f 1 ) ) 

check to see if contour line intersects edge 
between point 1 and point 3 

if ((cl-fl)*(cl-f3) .le.O. ) then 

contour line intersects edges 
between point 1 and point 4 and 
between point 1 and point 3 

ix2 ( i , l)=ofsetx+scalex*(xl+(x3-xl)*(cl~fl)/(f3“fl) ) 
iy2 ( i , 1 )=of sety+scaley* (yl+(y3-yl)*(cl-fl)/(f3-fl) ) 
ind( i , 1 ) = . true . 

must be a problem 

else 

set error flag 
f lag( i )=. true . 
endif 

contour line did not intersect edges 
between point 3 and point 4 and 
between point 1 and point 4 

check to see if contour line intersects edge 
between point 1 and point 3 

else if ( (cl-fl)*(cl-f3) . le.O. ) then 

contour line only intersects edge 
between point 1 and point 2 
must be a problem 
set error flag 

f lag( i )= . true . 

endif 

endif 

endif 
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500 continue 
c 

do 550 i=is f iel 
c 

c check to see if an error occured 

c 

if (flag(i)) then 
c 

c issue error message and stop 

c 

write(*,*) 'There is something fishy about conplt!' 
stop 
c 

endif 

c 

550 continue 
c 

do 560 i=is,iel 
c 

c check to see if line is to be plotted 

c 

if ( ind( i , 1 ) ) then 
c 

c plot a single line 

c 

call ufblin( ixl ( i , 1 ) t iyl( i t 1 ) , ix2 ( i , 1 ) , iy2 ( i f 1 ) ) 
c 

c check for more lines 

c 

if ( ind( i , 2 ) ) then 
c 

c plot two more lines 

c 

call ufblin(ixl(i, 2) , iyl(i,2) , ix2(i,2) , iy2(i,2) ) 
call ufblin(ixl(i, 3) , iyl ( i # 3 ) , 1x2 ( i , 3 ) , iy2(i,3) ) 
c 

endif 

c 

endif 

c 

560 continue 
c 

400 continue 
c 

c increment contour level 

c 

cl=cl+delf 

c 

c plot contours 

c 

go to 150 
c 

end 

c************************************* ******************************* 

subroutine fminmax( f , fmin, fmax .is.ie.js.je^m.n.iflag) 
c 



UO + 


33 


c 

c 

c 

c 

c 


c 

c 


c 

c 

c 


c 

100 

c 


dimension f(m,n) 
data eps/lel4/ 

check to reset minimums and maximums 

if (iflag.eq.O) then 

fmin=eps 

fmax=-eps 

end i f 

do 100 j=js,je 
do 100 i=is,ie 

check to reset minimums and maximums 

fmax=amaxl ( f ( i , j ) , fmax ) 
fmin=aminl ( f ( i , j ) , fmin) 

continue 

return 

end 




C 

subroutine ge tf ( td , tu , tv , te , dd , m , n , i var ) 
c 

dimension td(*), tu(*), tv(*), te(*), dd(*) 
c 

common/f luid/gamma, gml , gpl , gmlg, gplg, ggml 
c 

data cv/4390./ 
c 

mn=m*n 

c 

rcv=l./cv 

rggml=l . /ggml 

coel=(gamma-l . 0)*cv 

coe2=gamma* ( gamma- 1 . 0 ) *cv 

c 

if ( ivar . eq . 1 ) then 
do 10 i=l,mn 
dd( i )=td( i ) 

10 continue 

return 

end if 
c 

i f ( ivar . eq . 2 ) then 
do 20 i=l,mn 
rrho=l . /td( i ) 
ul=tu( i ) *rrho 
vl=tv( i ) *rrho 



dd( i ) = ( te( i ) *rrho-0 . 5* (ul*ul+vl*vl ) ) *rcv 
continue 
return 
end if 

if ( ivar . eq. 3 ) then 
do 30 i=l,mn 
rrho=l./td(i) 
ul=tu( i ) *rrho 
vl=tv( i ) *rrho 

tl=( te( i ) *rrho-0 . 5* (ul*ul+vl*vl ) )*rcv 
dd( i )=td( i ) *gml*cv* tl 
continue 
return 
end if 

if ( ivar. eq. 4) then 
do 40 i=l,mn 

rrho=l . /td( i ) 
ul=tu ( i ) *rrho 
vl=tv( i ) *rrho 
uv=ul*ul+vl*vl 
tl=( te ( i ) *rrho-0 . 5*uv) *rcv 
cl=abs ( tl*ggml*cv) 
dd( i )=sqrt (uv/cl ) 
continue 
end if 

if (ivar. eq. 5) then 
do 50 i=l,mn 
rrho=l./td(i) 
ul=tu ( i ) * rrho 
vl=tv( i ) *rrho 

tl=( te( i ) *rrho-0 . 5* (ul*ul+vl*vl ) )*rcv 
pl=td( i ) *gml*cv*tl 
dd( i )=log( tl**rggml/pl ) 
continue 
end if 

if ( ivar . eq. 6) then 
do 60 i=l,mn 

uv= tu ( i ) * tu ( i ) + tv ( i ) * tv ( i ) 
dd( i )=0 . 5*uv/td( i ) 
continue 
end if 

i f ( ivar . eq . 7 ) then 
do 70 i=l,mn 
rrho=l./td(i) 
ul=tu( i ) *rrho 
vl=tv( i ) *rrho 

tl=( te( i ) *rrho-0 .5*(ul*ul*vl*vl) )*rcv 
pl=td( i ) *gml*cv* tl 
dd( i )=( te( i ) +pl ) *rrho 
continue 
end if 
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c 

return 

end 

subroutine gridxy(x,y,m,n,is,ie,js,je,iflag) 
c 

dimension x(m,n) , y(m, n) 
c 

common /maping/ scalex, ofsetx, scaley, ofsety, scaled, ofsetd 
real scalex, ofsetx, scaley, ofsety, scaled, ofsetd 

c 

dimension ix(4),iy(4) 
c 

c iflag=0 draw the boundaries only 

c iflag=l draw constant-i lines only 

c iflag=2 draw constant- j lines only 

c iflag=3 draw the entire gird 

c 

c set grid line look up table (lut) index (color) 

c 

call ufbcol(255) 
c 

c compute do loop end and increment 

c 

iil=ie-l 
j jl=je-l 
ii=l 

JJ=1 

c 

c check to reset increments 

c 

if(iflag.eq.O.or.iflag.eq.l) j j=je-js 
if(iflag.eq.0.or.iflag.eq.2) ii=ie-is 
c 

c check to plot grid lines 

c 

if ( ( i f lag . ne . 2 ) . and . ( i i . gt . 0 ) . and . ( j s . ne . j e ) ) then 
c 

do 100 i=is,ie,ii 

do 100 j=js,jjl 

c 

c convert world coordinates to screen coordinates 

c 

ix( 1 )=scalex*x( i , j )+ofsetx 
ix(2 )=scalex*x( i,j+l)+ofsetx 
iy( 1 )=scaley*y( i , j )+ofsety 
iy(2 )=scaley*y( i,j+l)+ofBety 
c 

c plot grid line 

c 

call ufblin( ix( l),iy(l),ix(2),iy(2)) 
c 

100 continue 

c 

endif 

c 
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c check to plot grid lines 

c 

lf((iflag.ne.l). and . ( j j . gt . 0 ) . and . ( i s . ne . ie ) ) then 
c 

do 200 j=js,je,jj 

do 200 i=is, ill 

c 

c convert world coordinates to screen coordinates 

c 

ix( 1 )=scalex*x( i ,j)+ofsetx 
ix(2)=scalex*x(i+l,j) +of setx 
iy( 1 )=scaley*y( i ,j)+ofsety 
iy(2)=scaley*y(i+l,j)+ofsety 
c 

c plot grid line 

c 

call ufblin( ix(l),iy(l),ix(2),iy(2)) 
c 

200 continue 

c 

endif 

c 

return 

end 


c 

subroutine linmap( al , a2 , bl , b2 , cl , c2 ) 
c 

c compute scale factor and offset 

c 

if(a2.ne.al) then 

cl=(b2-bl)/(a2-al) 
c2=bl - cl*al + .5 

else 

cl=0. 

c2=(bl+b2 ) *0 . 5 
end if 
return 
end 

★ iHr ★ "k 


C 

C 


C 


C 


subroutine plotd( x , y, q, d , m, n , ivar , i s , 1 e , j s , je , nc , sides , 
* id, idmin, idmax , i t, i f 1 ag) 


dimension x(m,n),y(m,n),q(m,n,4),d(m.n) 


common /extrem/ xmin, xmax, ymin, 

real xmin, xmax, ymin. 


ymax, dmin, dmax 

ymax, dmin, dmax 


common /maping/ scalex, 
real scalex. 


ofsetx, scaley, 
ofsetx, scaley. 


ofsety, scaled, ofsetd 
ofsety, scaled, ofsetd 


c 


common/fluid/gamma, gml, gpl , gmlg, gplg, ggml 
real idmin, idmax 
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c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


save odmin, odmax 

data gamma, odmin, odmax/1 . 4, +le30, -le30/ 
contour plot 

written by: Eric L. Everton & Jamshid S. Abolhassani 

gml=gamma- 1 . 
gpl=gamma+ 1 . 
gmlg=gml/gamma 
gplg=gpl/gamma 
ggml=gamma*gml 

ivar=0 write out the grid 

1 density 

2 temperature 

3 pressure 

4 mach number 

5 entropy 

6 dynamic pressure 

7 total enthalpy 

scaling the grid & flow variables 
if (it.eq.O) then 

compute x and y minimums and maximums 

call fminmax ( x , xmin, xmax ,is,ie,js,je,m,n,0) 
call fminmax ( y, ymin, ymax, is, le,js,je,m,n,0) 

compute aspect ratio and 
x and y ranges 

yoverx=1024 . /1280 . 

arange=xmax-xmin 

brange=ymax-ymin 

check for zero ranges 

if ( ( arange . eq. 0. ) . or . (brange . eq. 0 . ) ) then 
a range is zero 

issue an error message and stop 

wr; ^ e ( * » * ) xmin - xmax or ymin - ymax ea O' 
stop 

endif 

check aspect ratio with range ratio 

if (brange/arange. gt. yoverx) then 

adjust xmin and xmax 
recompute range 
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c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


c 


hafdif- .5 * (brange/yoverx - arange) 
xmin =xmin-hafdif 
xmax =xmax+hafdi f 
arange=xmax-xmin 

else 

adjust ymin and ymax 
recompute range 

hafdif= .5 * ( arange*yoverx - brange) 
ymin =ymin-hafdif 
ymax =ymax+hafdif 
brange=ymax-ymin 

endif 

add border 

xmin=xmin - sides*arange 
xmax=xmax + sides*arange 
ymin=ymin - sides*brange 
ymax=ymax + sides*brange 

compute scale and offset for x and y 

ca ^ ^ 1 inmap ( xmin , xmax , 0 . , 1279. , sea lex , of set x ) 
call linmap( ymin, ymax, 0. , 1023 . , scaley, of sety ) 

endif 

get d variable 

call getf (q( 1 , 1 , 1 ) , q( i f 1 f 2 ) , q( 1 , 1 , 3 ) , q( i , 1 , 4) , d, m, n 

compute d minimum and maximum 

call fminmax(d, edmin, edmax, is,ie,js,je,m,n, 0 ) 
idflag=l 

if (edmin. eq. edmax) idflag=0 

save overall d minimum and maximum 

if ( edmin. It . odmin) odmin=cdmin 
if ( edmax . gt . odmax ) odmax=cdmax 

check to use current d minimum and maximum 

if (id.eq.O) then 

set dmin and dmax to current 

dmin=cdmin 

dmax=cdmax 


ivar ) 



else 
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c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


set dmln and dmax to Input 

dmin=idmin 

dmax=idmax 

endif 

compute scale and offset for d 
if ( idf lag. eq. 1 ) 

1 call 1 inmap ( dmin , dmax , 2 . ,255. f scaled , of setd ) 

write current, overall & used d minimums and maximums 

write(62,*) ' current dmin & dmax = cdmin, cdmax 

write(62,*) ' overall dmin & dmax = odmin, odmax 

write ( 62 , * ) ' used dmin & dmax = dmin, dmax 

clear frame buffer 

call ufbcle 

write frame number 

write(62,*)' Plotting frame number ' , i t 
plot contours 
if (idflag.eq. 1) 

1 call conplt(x,y,d,m,n, is, ie, js, je,2. ,255. ,nc) 
plot grid lines 

call gridxy (x, y , m , n, is , ie , j s , j e , if lag) 

send frame buffer to display 

call ufbwri 

return 

end 




